%% Calculate pi_primes: exports, intermediate inputs, labor

% function [pi_prime_nkjf,pi_prime_nkj, pi_cprime_inj] =...
%     TradeShareCES_fun(P_hat_inj,gamma,w_hat,pi_inj,...
%     GAMMA_sorted,alpha_f_sorted,pi_nkjp_sorted,pi_c_inj)

function [pi_mnj1,pi_mnjf1,pi_c_mnj1,pi_M1,pi_l1,pi_M_f1,pi_l_f1] = ...
    TradeShareCESOlig_fun(P_hat_mnj,b_hat,b_hat_f,P_M_f_hat,pi_M0,pi_l0,w_hat,pi_mnj0,...
    pi_M_f0,pi_l_f0,pi_mnjf0,pi_c_mnj0,mu_hat_mnjf)
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Produced for JdG,AL,IM by Christopher Evans at UPF 
%
% Solves the new trade shares that are produced from taking into account 
% the different unit cost of firms, prices and wages. Trade shares are 
% needed in the expenditure X equation.
%
% NEW: For all code, we will change the subscripts to match those in
% master_notes3.pdf and the draft. I.e., we will change to {mn,ij} pairs,
% which represent country pair {mn} and sector pair {ij}. Global variables 
% changed and slight modification to file structure for calling up data.
% We also get rid of gamma_type variable. What previously _2.csv. files
% slightly as well as adjusting in our new directory.
%
% MAJOR CHANGE: Solution of price index now based on CES cost function,
% which implies a change of in both intermedite input and labor shares 
% given a shock. Therefore, the equilibrium outcomes will increase and more
% output from the function relative to before. 
% Will refer to master_notes3 for formula and solve for.
%
% This function will also be used for both the baseline setting as well as
% the heterogeneous calibration. Therefore, the initial labor and
% intermediate import shares will have to adjust away from the data (alpha,
% gamma) in the inital wedge exercise. Therefore, the function is set up to
% run with a pi_l and pi_M rather than data shares.
%
% We also get rid of using the .(ISO{}) cell call up for the gammas and
% instead use a large matrix saved in Step2. This will help save on loops
% and speed up the proceedure.
%
% Last Updated: 5/02/2019 by IM & JDG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global alphaT rhoT lambdaT etaT rho sigma lambda eta varphi ...
    tol maxit Pi_hat_n D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    N J FI_J france countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO ... 
    sL_n sPi_n sD_n alpha_WIOD_j alpha_WIOD_francej labour_share_PWT
    
% Initialize outcome variables to initial valued fed into function

pi_mnj = pi_mnj0;
pi_mnjf = pi_mnjf0;
pi_c_mnj = pi_c_mnj0;

pi_M = pi_M0;
pi_l = pi_l0;
pi_M_f = pi_M_f0;
pi_l_f = pi_l_f0;

%% Export shares

% pi_mnj1:  Export share by firms for each country given a representative
%           firm in a sector
%
%           Since we do not have firm level data where firm export
%           share must sum to 1 over firms for the country-sector 
%           exports we cannot update the sector level export share
%           which are equal to 1 for tradable sector and 0 for
%           non-tradable sectors except NT-NT trade within the same
%           country. Hence the time t+1 variable is equal to its time t
%           equivalent. 

pi_mnj1 = pi_mnj; 

%  pi_mnjf1 : Export share data for french firms

% We know for n=france, hence we construct the change in french firm export
% share. Based on master_notes3.pdf equations

rho_f=firmsecD_sorted*rho;

%cost_f = (mu_hat_mnjf.*repmat(b_hat_f,[1 N]).*repmat(a_hat_f,[1 N])).^(1-repmat(rho_f,[1 N]));

NUM_f = pi_mnjf.*Xi_hat_mnjf.*((mu_hat_mnjf.*repmat(b_hat_f,[1 N]).*repmat(a_hat_f,[1 N])).^(1-repmat(rho_f,[1 N])));

% For denominator use destination-specific price index raised to 
% (1-rho)

DEN_f = firmsecD_sorted*P_hat_mnj(start(france):start(france+1)-1,:).^(1-repmat(rho,[1,N]));

pi_mnjf1= NUM_f./DEN_f;
% disp('pi_mnjf1')
% [min(min(pi_mnjf1)) max(max(pi_mnjf1))] 
clear NUM_f DEN_f

% pi_c_mnj1: update country shares

% Need to expand out sigma to be able to multiply easily

sigma_nj_repmat=repmat(sigma,[N 1]);
sigma_mnj_repmat=repmat(sigma_nj_repmat,[1 N]);

% P_hat_pi_c: Calculating the numerator for the pi_c_mnj1 equation                      

%P_hat_pi_c = P_hat_mnj.^(1-sigma_mnj_repmat).*pi_c_mnj;

% sum_P_hat_pi_c: Summing to create the denominator of the pi__c_mnj 
%                 equation

sum_P_hat_pi_c= sum_by_country_dummy*(P_hat_mnj.^(1-sigma_mnj_repmat).*pi_c_mnj) ;

% We need to make sure there are no zeros because we will be dividing by
% this matrix replicated to N countries

sum_P_hat_pi_c=sum_P_hat_pi_c+0.000000001.*(sum_P_hat_pi_c==0);

% sum_P_hat_pi_c_repmat: replicate the matrix for each i country

%sum_P_hat_pi_c_repmat=repmat(sum_P_hat_pi_c,[N 1]);

% pi_c_mnj1: Share of final consumption spending on goods from m in the
%          total consumption spending on goods in sector j, country,
%           n at time period t+1 (after the shock)

pi_c_mnj1 = (P_hat_mnj.^(1-sigma_mnj_repmat).*pi_c_mnj)./(repmat(sum_P_hat_pi_c,[N 1]));
% disp('pi_c_mnj1')
% [min(min(pi_c_mnj1)) max(max(pi_c_mnj1))] 
%% Labor shares

% pi_l0: solve for value based on data labor shares and updated wage and
%        price vectors


% pi_l1: labor share for representative firms across countries

% Taking wage of france to power of 1-lambda and multiplying by labor share

%w_hat_to_vector=countrysecD*w_hat';
%NUM = pi_l.*w_hat_to_vector.^(1-lambda);

% Using b_hat for denominator

%DEN = (b_hat).^(1-lambda);

%pi_l1 = NUM./repmat(DEN,[1,N]);
pi_l1 = (pi_l.*(countrysecD*w_hat').^(1-lambda))./((b_hat).^(1-lambda));
% Get single vector of repeated matrix

%pi_l1 = pi_l1(:,1);

% disp('pi_l1')
% [min(pi_l1)' max(pi_l1)']

% Replace NAN with zeros 

pi_l1(isnan(pi_l1))=0;

% pi_l_f1: labor share at firm level, which is a function of previous share
%          as well as cost function and price level
    
% Taking wage of france to power of 1-lambda and multiplying by labor share
  
%NUM_f = pi_l_f.*w_hat(france).^(1-lambda);

%NUM_f = pi_l_f.*repmat(w_hat(france),[FI_J 1]).^(1-lambda);

% Using b_hat_f for denominator

%DEN_f = (b_hat_f).^(1-lambda);

pi_l_f1 = (pi_l_f.*repmat(w_hat(france),[FI_J 1]).^(1-lambda))./((b_hat_f).^(1-lambda));
% disp('pi_l_f1')
% [min(pi_l_f1)' max(pi_l_f1)']

%clear NUM DEN NUM_f DEN_f;

%% Import shares

% pi_M1: import share for representative firm across countries, which is a
%        function of pervious sahre as well as updated price indices that
%        have been fed in; i.e,. P_hat_mnj

% P_hat_mnj: extending the matrix by J for  prices
   
P_hat_mnj_repmat = kron(P_hat_mnj,ones(1,J));
   
% P_hat_weta = pi_M.*(P_hat_mnj_repmat)^(1-eta);

NUM = pi_M.*P_hat_mnj_repmat.^(1-eta);

% Summing by n and by i, so we just need to sum along the rows, this
% will give us a transposed vector [ 1xJ*N ] and the price of
% intermediate goods ==> P_M_hat except not raised to 1/(1-eta)
   
%DEN = squeeze(sum(NUM,1));

% Stretch out to divide by numerator, so repeat transpose across J*N

%DEN1 = repmat(DEN', [1 N*J]);
DEN1 = repmat(squeeze(sum(NUM,1)),[N*J 1]);

% Add very small number to insure not dividing by zero. Can occur given
% initial import shares of 0 being summed over

DEN1 = DEN1+.0000000001.*(DEN1==0);

pi_M1 = NUM./DEN1;

% disp('pi_M1')
% [min(min(pi_M1)) max(max(pi_M1))]

clear NUM DEN1;

% pi_M_f1: import share at firm level, which is a function of previous
%          share as well as update price indices that have been fed in;
%          i.e. P_hat_mnj

%DEN_f_repmat = repmat(P_M_f_hat.^(1-eta),[1,J]);
DEN_f_repmat = repmat(P_M_f_hat.^(1-eta),[1,N*J]);

% To create share we need to loop over country sources to generate the
% numerator and  divide by denominator. 

% for n = 1:N
%     P_hat_mnjf_0 = ...
%       repmat(P_hat_mnj(start(n):start(n+1)-1,france)',[FI_J 1]); 
%     pi_M_f1.(ISO{1,n}) = (pi_M_f.(ISO{1,n}).*P_hat_mnjf_0.^(1-eta))./DEN_f_repmat;
% end

%P_hat_mnjf_0 = repmat(P_hat_mnj(:,france)',[FI_J,1]);
pi_M_f1 = (pi_M_f.*(repmat(P_hat_mnj(:,france)',[FI_J,1])).^(1-eta))./DEN_f_repmat;

end